Introdução

“Infecção hospitalar é uma infecção adquirida após a admissão do paciente na unidade hospitalar e pode se manifestar durante a internação ou após a alta. Pela sua gravidade e aumento do tempo de internação do paciente, é causa importante de morbidade e mortalidade, caracterizando-se como problema de saúde pública.”

Os dados disponíveis no arquivo data/SCENIC.txt foram coletados pelo CDC (US Center for Disease Control), no âmbito do Projeto SCENIC. O principal objetivo do projeto era determinar se programas de vigilância e controle foram capazes de reduzir as taxas de infecção hospitalar. Os dados referem-se a uma amostra de 113 hospitais selecionados a partir de um conjunto de 338 hospitais avaliados. Cada linha do conjunto de dados contém uma identificação (1-113) e fornece informação a respeito de 11 variáveis para um único hospital. Os dados apresentados referem-se ao período de 1975-1976.

As 12 variáveis são:

  1. IDnumber: 1-113 (identificação do hospital)
  2. LengthStay: período de internação médio de todos os pacientes no hospital (em dias)
  3. Age: idade média dos pacientes (anos)
  4. InfectRisk: risco de infecção, calculado como a probabilidade média estimada de contrair infecção no hospital (em %)
  5. CultRatio: razão do número de culturas realizadas pelo número de pacientes sem sintomas de infecção hospitalar, vezes 100
  6. XrayRatio: razão do número de raios-X realizados pelo número de pacientes sem sintomas de pneumonia, vezes 100
  7. NBeds: número médio de leitos do hospital, durante o período avaliado
  8. MedSchool: afiliação a alguma Escola de Medicina (1=Sim, 2=Não)
  9. Region: região geográfica (1=NE, 2=NC, 3=S, 4=W)
  10. DailyCensus: número médio de pacientes no hospital por dia, durante o período avaliado
  11. NNurses: número médio de enfermeiros no hospital
  12. Facilities: percentual de 35 serviços providos pelo hospital

Acredita-se que o período médio de internação de um paciente LengthStay (variável de resposta) possa ser previsto a partir do risco de infecção hospitalar, bem como outras características do hospital e de procedimentos de rotina realizados.

Análise Exploratória de Dados

Conduza a análise exploratória da massa de dados SCENIC, a fim de compreender suas características principais.

rm(list=ls())
scenic <- read.table("data/SCENIC.txt", header = FALSE)

str(scenic)
## 'data.frame':    113 obs. of  12 variables:
##  $ V1 : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ V2 : num  7.13 8.82 8.34 8.95 11.2 ...
##  $ V3 : num  55.7 58.2 56.9 53.7 56.5 50.9 57.8 45.7 48.2 56.3 ...
##  $ V4 : num  4.1 1.6 2.7 5.6 5.7 5.1 4.6 5.4 4.3 6.3 ...
##  $ V5 : num  9 3.8 8.1 18.9 34.5 21.9 16.7 60.5 24.4 29.6 ...
##  $ V6 : num  39.6 51.7 74 122.8 88.9 ...
##  $ V7 : int  279 80 107 147 180 150 186 640 182 85 ...
##  $ V8 : int  2 2 2 2 2 2 2 1 2 2 ...
##  $ V9 : int  4 2 3 4 1 2 3 2 3 1 ...
##  $ V10: int  207 51 82 53 134 147 151 399 130 59 ...
##  $ V11: int  241 52 54 148 151 106 129 360 118 66 ...
##  $ V12: num  60 40 20 40 40 40 40 60 40 40 ...
colnames(scenic) <- c("IDnumber","LengthStay","Age","InfectRisk","CultRatio","XrayRatio","NBeds","MedSchool","Region","DailyCensus","NNurses","Facilities")

scenic$MedSchool[scenic$MedSchool == 2 ] <- 0

names(scenic)
##  [1] "IDnumber"    "LengthStay"  "Age"         "InfectRisk"  "CultRatio"  
##  [6] "XrayRatio"   "NBeds"       "MedSchool"   "Region"      "DailyCensus"
## [11] "NNurses"     "Facilities"
summary(scenic)
##     IDnumber     LengthStay          Age          InfectRisk      CultRatio    
##  Min.   :  1   Min.   : 6.700   Min.   :38.80   Min.   :1.300   Min.   : 1.60  
##  1st Qu.: 29   1st Qu.: 8.340   1st Qu.:50.90   1st Qu.:3.700   1st Qu.: 8.40  
##  Median : 57   Median : 9.420   Median :53.20   Median :4.400   Median :14.10  
##  Mean   : 57   Mean   : 9.648   Mean   :53.23   Mean   :4.355   Mean   :15.79  
##  3rd Qu.: 85   3rd Qu.:10.470   3rd Qu.:56.20   3rd Qu.:5.200   3rd Qu.:20.30  
##  Max.   :113   Max.   :19.560   Max.   :65.90   Max.   :7.800   Max.   :60.50  
##    XrayRatio          NBeds         MedSchool          Region     
##  Min.   : 39.60   Min.   : 29.0   Min.   :0.0000   Min.   :1.000  
##  1st Qu.: 69.50   1st Qu.:106.0   1st Qu.:0.0000   1st Qu.:2.000  
##  Median : 82.30   Median :186.0   Median :0.0000   Median :2.000  
##  Mean   : 81.63   Mean   :252.2   Mean   :0.1504   Mean   :2.363  
##  3rd Qu.: 94.10   3rd Qu.:312.0   3rd Qu.:0.0000   3rd Qu.:3.000  
##  Max.   :133.50   Max.   :835.0   Max.   :1.0000   Max.   :4.000  
##   DailyCensus       NNurses        Facilities   
##  Min.   : 20.0   Min.   : 14.0   Min.   : 5.70  
##  1st Qu.: 68.0   1st Qu.: 66.0   1st Qu.:31.40  
##  Median :143.0   Median :132.0   Median :42.90  
##  Mean   :191.4   Mean   :173.2   Mean   :43.16  
##  3rd Qu.:252.0   3rd Qu.:218.0   3rd Qu.:54.30  
##  Max.   :791.0   Max.   :656.0   Max.   :80.00

Transformação do tipo de variáveis

scenic$IDnumber <- as.factor(scenic$IDnumber)
scenic$MedSchool <- as.factor(scenic$MedSchool)
scenic$Region <- as.factor(scenic$Region)

Para utilizar os nomes das variáveis diretamente

attach(scenic)

Boxplot e histograma com linha de densidade estimada da variável de resposta.

boxplot(LengthStay, horizontal = TRUE)

hist(LengthStay, freq = FALSE, breaks = "FD", main = "")
lines(density(LengthStay), col="blue", lwd=3)

#matriz de gráficos de dispersão
plot(scenic[, -c(1,8,9)])

#carregar biblioteca para fazer matriz de correlação
library(corrplot)
## corrplot 0.84 loaded
res <- cor(scenic[, -c(1,8,9)])
round(res,2)
##             LengthStay   Age InfectRisk CultRatio XrayRatio NBeds DailyCensus
## LengthStay        1.00  0.19       0.53      0.33      0.38  0.41        0.47
## Age               0.19  1.00       0.00     -0.23     -0.02 -0.06       -0.05
## InfectRisk        0.53  0.00       1.00      0.56      0.45  0.36        0.38
## CultRatio         0.33 -0.23       0.56      1.00      0.42  0.14        0.14
## XrayRatio         0.38 -0.02       0.45      0.42      1.00  0.05        0.06
## NBeds             0.41 -0.06       0.36      0.14      0.05  1.00        0.98
## DailyCensus       0.47 -0.05       0.38      0.14      0.06  0.98        1.00
## NNurses           0.34 -0.08       0.39      0.20      0.08  0.92        0.91
## Facilities        0.36 -0.04       0.41      0.19      0.11  0.79        0.78
##             NNurses Facilities
## LengthStay     0.34       0.36
## Age           -0.08      -0.04
## InfectRisk     0.39       0.41
## CultRatio      0.20       0.19
## XrayRatio      0.08       0.11
## NBeds          0.92       0.79
## DailyCensus    0.91       0.78
## NNurses        1.00       0.78
## Facilities     0.78       1.00
corrplot(res, method = "circle", order = "hclust", sig.level = 0.01, insig = "blank")

A matriz obtida sugere correlação entre algumas variáveis, cujos diagramas de dispersão serão exibidos a seguir. Além disso, Scatterplot para variável de resposta em função das explicativas filtrado por categóricas.

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
  #MedSchool
scatterplot(LengthStay ~ DailyCensus|MedSchool, data = scenic)

scatterplot(LengthStay ~ InfectRisk|MedSchool, data = scenic)

scatterplot(LengthStay ~ Facilities|MedSchool, data = scenic)

  #Region
scatterplot(LengthStay ~ CultRatio|Region, data = scenic)

scatterplot(LengthStay ~ XrayRatio|Region, data = scenic)

scatterplot(LengthStay ~ Facilities|Region, data = scenic)

scatterplot(LengthStay ~ InfectRisk|Region, data = scenic)

#Construção de nova base de dados
detach(scenic)
scenic <- scenic[, -c(1)]
attach(scenic)
data_rls <- data.frame(LengthStay = LengthStay, InfectRisk = InfectRisk, Facilities = Facilities, XrayRatio = XrayRatio, Age = Age, CultRatio =  CultRatio, NBeds =  NBeds, MedSchool = MedSchool, Region = Region, DailyCensus = DailyCensus, NNurses = NNurses)

summary(data_rls)
##    LengthStay       InfectRisk      Facilities      XrayRatio     
##  Min.   : 6.700   Min.   :1.300   Min.   : 5.70   Min.   : 39.60  
##  1st Qu.: 8.340   1st Qu.:3.700   1st Qu.:31.40   1st Qu.: 69.50  
##  Median : 9.420   Median :4.400   Median :42.90   Median : 82.30  
##  Mean   : 9.648   Mean   :4.355   Mean   :43.16   Mean   : 81.63  
##  3rd Qu.:10.470   3rd Qu.:5.200   3rd Qu.:54.30   3rd Qu.: 94.10  
##  Max.   :19.560   Max.   :7.800   Max.   :80.00   Max.   :133.50  
##       Age          CultRatio         NBeds       MedSchool Region
##  Min.   :38.80   Min.   : 1.60   Min.   : 29.0   0:96      1:28  
##  1st Qu.:50.90   1st Qu.: 8.40   1st Qu.:106.0   1:17      2:32  
##  Median :53.20   Median :14.10   Median :186.0             3:37  
##  Mean   :53.23   Mean   :15.79   Mean   :252.2             4:16  
##  3rd Qu.:56.20   3rd Qu.:20.30   3rd Qu.:312.0                   
##  Max.   :65.90   Max.   :60.50   Max.   :835.0                   
##   DailyCensus       NNurses     
##  Min.   : 20.0   Min.   : 14.0  
##  1st Qu.: 68.0   1st Qu.: 66.0  
##  Median :143.0   Median :132.0  
##  Mean   :191.4   Mean   :173.2  
##  3rd Qu.:252.0   3rd Qu.:218.0  
##  Max.   :791.0   Max.   :656.0

Note que as variáveis LengthStay está relacionada com InfectRisk e XrayRatio, contudo não linear. Portanto vamos analisar os plots de maneira separada para decidir se cabe a criação de um modelo de ordem superior para analisar a significância da variável na resposta.

plot(LengthStay ~ InfectRisk)

plot(LengthStay ~ XrayRatio)

testIR_lm <- lm(LengthStay ~ I(InfectRisk*InfectRisk), data = data_rls)
testXR_lm <- lm(LengthStay ~ I(XrayRatio*XrayRatio), data = data_rls)

summary(testIR_lm)
## 
## Call:
## lm(formula = LengthStay ~ I(InfectRisk * InfectRisk), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9031 -0.8445 -0.0828  0.5875  7.9534 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 7.75885    0.30447  25.483  < 2e-16 ***
## I(InfectRisk * InfectRisk)  0.09107    0.01278   7.125  1.1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.59 on 111 degrees of freedom
## Multiple R-squared:  0.3139, Adjusted R-squared:  0.3077 
## F-statistic: 50.77 on 1 and 111 DF,  p-value: 1.104e-10
summary(testXR_lm)
## 
## Call:
## lm(formula = LengthStay ~ I(XrayRatio * XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8791 -1.1153 -0.2815  0.8144  8.4927 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              7.954e+00  4.005e-01  19.863  < 2e-16 ***
## I(XrayRatio * XrayRatio) 2.408e-04  5.185e-05   4.644 9.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.757 on 111 degrees of freedom
## Multiple R-squared:  0.1627, Adjusted R-squared:  0.1552 
## F-statistic: 21.57 on 1 and 111 DF,  p-value: 9.435e-06

Nota-se que apesar dos coeficientes angulares dos modelos e os valores da estatística t serem significativos em ambos os casos, o valor de Raj^2 não foi tão expressivo. Contudo seguiremos com mais variáveis e deixaremos uma possível exclusão no método de Eliminação Regressiva.

Vamos iniciar utilizando o método de Eliminação Regressiva (Backward Elimination). O procedimento começa com o ajuste do modelo completo e, a cada etapa, eliminamos a variável explicativa com maior valor-p acima do nível de significância 0.01:

# modelo completo
m.all <- lm(LengthStay ~. +I(InfectRisk*InfectRisk)+I(XrayRatio*XrayRatio)  , data = data_rls)
summary(m.all)
## 
## Call:
## lm(formula = LengthStay ~ . + I(InfectRisk * InfectRisk) + I(XrayRatio * 
##     XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2512 -0.7058 -0.0711  0.5724  5.9347 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 6.5566865  2.4280121   2.700 0.008159 ** 
## InfectRisk                 -0.2284972  0.4932714  -0.463 0.644228    
## Facilities                 -0.0070984  0.0139253  -0.510 0.611371    
## XrayRatio                  -0.0303683  0.0452581  -0.671 0.503798    
## Age                         0.0654384  0.0288074   2.272 0.025298 *  
## CultRatio                   0.0041990  0.0161557   0.260 0.795478    
## NBeds                      -0.0058389  0.0035851  -1.629 0.106592    
## MedSchool1                  0.1138183  0.4423737   0.257 0.797494    
## Region2                    -0.7676628  0.3494305  -2.197 0.030385 *  
## Region3                    -1.0931318  0.3546769  -3.082 0.002670 ** 
## Region4                    -1.7531293  0.4436883  -3.951 0.000147 ***
## DailyCensus                 0.0161506  0.0043833   3.685 0.000375 ***
## NNurses                    -0.0056739  0.0022501  -2.522 0.013294 *  
## I(InfectRisk * InfectRisk)  0.0760331  0.0544231   1.397 0.165547    
## I(XrayRatio * XrayRatio)    0.0002723  0.0002715   1.003 0.318212    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.213 on 98 degrees of freedom
## Multiple R-squared:  0.6474, Adjusted R-squared:  0.5971 
## F-statistic: 12.85 on 14 and 98 DF,  p-value: < 2.2e-16

Vamos iniciar pela variável MedSchool, que apresenta maior valor-p (0.797), de acordo com o resumo da regressão para o modelo completo.

m.be <- update(m.all, . ~. - MedSchool)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ InfectRisk + Facilities + XrayRatio + 
##     Age + CultRatio + NBeds + Region + DailyCensus + NNurses + 
##     I(InfectRisk * InfectRisk) + I(XrayRatio * XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2448 -0.7188 -0.0625  0.5654  5.9096 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 6.5938253  2.4122598   2.733 0.007425 ** 
## InfectRisk                 -0.2558594  0.4793938  -0.534 0.594735    
## Facilities                 -0.0065794  0.0137133  -0.480 0.632438    
## XrayRatio                  -0.0301223  0.0450341  -0.669 0.505130    
## Age                         0.0648303  0.0285746   2.269 0.025450 *  
## CultRatio                   0.0050549  0.0157348   0.321 0.748695    
## NBeds                      -0.0060224  0.0034969  -1.722 0.088150 .  
## Region2                    -0.7570943  0.3453674  -2.192 0.030713 *  
## Region3                    -1.0926382  0.3529951  -3.095 0.002557 ** 
## Region4                    -1.7307990  0.4330600  -3.997 0.000124 ***
## DailyCensus                 0.0164860  0.0041652   3.958 0.000142 ***
## NNurses                    -0.0056280  0.0022324  -2.521 0.013298 *  
## I(InfectRisk * InfectRisk)  0.0784556  0.0533492   1.471 0.144569    
## I(XrayRatio * XrayRatio)    0.0002714  0.0002701   1.004 0.317605    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.208 on 99 degrees of freedom
## Multiple R-squared:  0.6472, Adjusted R-squared:  0.6009 
## F-statistic: 13.97 on 13 and 99 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é CultRatio
m.be <- update(m.be, . ~. - CultRatio)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ InfectRisk + Facilities + XrayRatio + 
##     Age + NBeds + Region + DailyCensus + NNurses + I(InfectRisk * 
##     InfectRisk) + I(XrayRatio * XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2367 -0.7721 -0.0447  0.5721  5.8764 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 6.7147998  2.3719782   2.831 0.005612 ** 
## InfectRisk                 -0.2730700  0.4742501  -0.576 0.566049    
## Facilities                 -0.0067007  0.0136465  -0.491 0.624485    
## XrayRatio                  -0.0270895  0.0438355  -0.618 0.537994    
## Age                         0.0617717  0.0268208   2.303 0.023342 *  
## NBeds                      -0.0059125  0.0034644  -1.707 0.090995 .  
## Region2                    -0.7834009  0.3340114  -2.345 0.020980 *  
## Region3                    -1.1194756  0.3414262  -3.279 0.001435 ** 
## Region4                    -1.7647991  0.4180412  -4.222 5.35e-05 ***
## DailyCensus                 0.0162553  0.0040844   3.980 0.000131 ***
## NNurses                    -0.0055141  0.0021942  -2.513 0.013567 *  
## I(InfectRisk * InfectRisk)  0.0825113  0.0516009   1.599 0.112969    
## I(XrayRatio * XrayRatio)    0.0002554  0.0002644   0.966 0.336290    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.202 on 100 degrees of freedom
## Multiple R-squared:  0.6468, Adjusted R-squared:  0.6045 
## F-statistic: 15.26 on 12 and 100 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é Facilities
m.be <- update(m.be, . ~. - Facilities)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ InfectRisk + XrayRatio + Age + NBeds + 
##     Region + DailyCensus + NNurses + I(InfectRisk * InfectRisk) + 
##     I(XrayRatio * XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2519 -0.7391 -0.0750  0.5324  5.8555 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 6.6075739  2.3530146   2.808  0.00598 ** 
## InfectRisk                 -0.3285045  0.4588829  -0.716  0.47572    
## XrayRatio                  -0.0256150  0.0435679  -0.588  0.55789    
## Age                         0.0610403  0.0266786   2.288  0.02422 *  
## NBeds                      -0.0063573  0.0033314  -1.908  0.05919 .  
## Region2                    -0.7704666  0.3317178  -2.323  0.02220 *  
## Region3                    -1.0826953  0.3318543  -3.263  0.00151 ** 
## Region4                    -1.7300329  0.4104508  -4.215 5.45e-05 ***
## DailyCensus                 0.0164760  0.0040443   4.074 9.21e-05 ***
## NNurses                    -0.0056765  0.0021609  -2.627  0.00996 ** 
## I(InfectRisk * InfectRisk)  0.0878634  0.0502467   1.749  0.08339 .  
## I(XrayRatio * XrayRatio)    0.0002475  0.0002629   0.942  0.34867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.198 on 101 degrees of freedom
## Multiple R-squared:  0.646,  Adjusted R-squared:  0.6074 
## F-statistic: 16.75 on 11 and 101 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é XrayRatio
m.be <- update(m.be, . ~. - XrayRatio)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ InfectRisk + Age + NBeds + Region + 
##     DailyCensus + NNurses + I(InfectRisk * InfectRisk) + I(XrayRatio * 
##     XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3098 -0.7744 -0.0726  0.5271  5.9073 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 5.742e+00  1.830e+00   3.138 0.002225 ** 
## InfectRisk                 -4.369e-01  4.189e-01  -1.043 0.299468    
## Age                         6.227e-02  2.651e-02   2.349 0.020761 *  
## NBeds                      -6.338e-03  3.320e-03  -1.909 0.059112 .  
## Region2                    -7.862e-01  3.296e-01  -2.385 0.018907 *  
## Region3                    -1.113e+00  3.267e-01  -3.408 0.000938 ***
## Region4                    -1.707e+00  4.072e-01  -4.191 5.92e-05 ***
## DailyCensus                 1.641e-02  4.030e-03   4.072 9.21e-05 ***
## NNurses                    -5.456e-03  2.121e-03  -2.572 0.011552 *  
## I(InfectRisk * InfectRisk)  9.920e-02  4.625e-02   2.145 0.034355 *  
## I(XrayRatio * XrayRatio)    9.490e-05  4.136e-05   2.295 0.023810 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.194 on 102 degrees of freedom
## Multiple R-squared:  0.6448, Adjusted R-squared:  0.6099 
## F-statistic: 18.51 on 10 and 102 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é InfectRisk
m.be <- update(m.be, . ~. - InfectRisk)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ Age + NBeds + Region + DailyCensus + 
##     NNurses + I(InfectRisk * InfectRisk) + I(XrayRatio * XrayRatio), 
##     data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3187 -0.7407 -0.0855  0.5983  5.9795 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4.626e+00  1.485e+00   3.115 0.002383 ** 
## Age                         6.720e-02  2.610e-02   2.575 0.011443 *  
## NBeds                      -6.022e-03  3.308e-03  -1.820 0.071620 .  
## Region2                    -7.798e-01  3.297e-01  -2.365 0.019885 *  
## Region3                    -1.091e+00  3.261e-01  -3.346 0.001147 ** 
## Region4                    -1.770e+00  4.029e-01  -4.393 2.72e-05 ***
## DailyCensus                 1.611e-02  4.021e-03   4.006 0.000117 ***
## NNurses                    -5.836e-03  2.091e-03  -2.791 0.006262 ** 
## I(InfectRisk * InfectRisk)  5.249e-02  1.158e-02   4.534 1.57e-05 ***
## I(XrayRatio * XrayRatio)    8.971e-05  4.107e-05   2.184 0.031234 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.194 on 103 degrees of freedom
## Multiple R-squared:  0.641,  Adjusted R-squared:  0.6096 
## F-statistic: 20.43 on 9 and 103 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é NBeds
m.be <- update(m.be, . ~. - NBeds)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ Age + Region + DailyCensus + NNurses + 
##     I(InfectRisk * InfectRisk) + I(XrayRatio * XrayRatio), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4170 -0.7893 -0.0832  0.5882  6.0746 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4.743e+00  1.500e+00   3.161 0.002057 ** 
## Age                         6.505e-02  2.636e-02   2.468 0.015218 *  
## Region2                    -9.176e-01  3.244e-01  -2.829 0.005609 ** 
## Region3                    -1.199e+00  3.242e-01  -3.700 0.000347 ***
## Region4                    -1.931e+00  3.973e-01  -4.861 4.15e-06 ***
## DailyCensus                 9.587e-03  1.846e-03   5.194 1.03e-06 ***
## NNurses                    -7.006e-03  2.012e-03  -3.483 0.000728 ***
## I(InfectRisk * InfectRisk)  5.383e-02  1.168e-02   4.608 1.16e-05 ***
## I(XrayRatio * XrayRatio)    8.966e-05  4.153e-05   2.159 0.033155 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.208 on 104 degrees of freedom
## Multiple R-squared:  0.6294, Adjusted R-squared:  0.6009 
## F-statistic: 22.08 on 8 and 104 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é XrayRatio*XrayRatio
m.be <- update(m.be, . ~. - I(XrayRatio*XrayRatio))
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ Age + Region + DailyCensus + NNurses + 
##     I(InfectRisk * InfectRisk), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7512 -0.8213 -0.1044  0.6629  6.3148 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 5.398374   1.494553   3.612 0.000468 ***
## Age                         0.063231   0.026801   2.359 0.020157 *  
## Region2                    -1.007715   0.327258  -3.079 0.002648 ** 
## Region3                    -1.378042   0.318864  -4.322 3.53e-05 ***
## Region4                    -2.101532   0.396166  -5.305 6.31e-07 ***
## DailyCensus                 0.009291   0.001873   4.961 2.71e-06 ***
## NNurses                    -0.006844   0.002045  -3.347 0.001136 ** 
## I(InfectRisk * InfectRisk)  0.063906   0.010895   5.866 5.25e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.228 on 105 degrees of freedom
## Multiple R-squared:  0.6128, Adjusted R-squared:  0.587 
## F-statistic: 23.74 on 7 and 105 DF,  p-value: < 2.2e-16
#A próxima variável a ser eliminada é Age
m.be <- update(m.be, . ~. - Age)
summary(m.be)
## 
## Call:
## lm(formula = LengthStay ~ Region + DailyCensus + NNurses + I(InfectRisk * 
##     InfectRisk), data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6662 -0.7627 -0.1206  0.7127  6.6221 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 8.823447   0.362821  24.319  < 2e-16 ***
## Region2                    -1.157607   0.327874  -3.531 0.000615 ***
## Region3                    -1.414600   0.325275  -4.349 3.15e-05 ***
## Region4                    -2.149807   0.404069  -5.320 5.81e-07 ***
## DailyCensus                 0.009504   0.001910   4.975 2.54e-06 ***
## NNurses                    -0.007263   0.002081  -3.490 0.000704 ***
## I(InfectRisk * InfectRisk)  0.065542   0.011104   5.902 4.36e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.255 on 106 degrees of freedom
## Multiple R-squared:  0.5923, Adjusted R-squared:  0.5692 
## F-statistic: 25.67 on 6 and 106 DF,  p-value: < 2.2e-16

Nesse ponto, todas as variáveis possuem valor-p menor que 0.01, então paramos com a iteração.Note que o coeficiente de determinação ajustado sofreu apenas uma pequena redução em comparação com aquele obtido para o modelo completo. Para o método de Eliminação Regressiva, as variáveis que devem ser usadas no modelo para LengthStay são Region, DailyCensus, NNurses e InfectRisk * InfectRisk e pelo Princípio da Hierarquia, como a variável InfectRisk * InfectRisk se manteve no modelo, devemos retornar com a variável InfectRisk.

Vamos proceder agora com o Método Gradativo de Seleção Sequencial (Stepwise Regression), que combina os procedimentos de Eliminação Regressiva e Seleção Progressiva:

# Stepwise Regression
library(leaps)
m.sr <- regsubsets(LengthStay ~. +I(InfectRisk*InfectRisk)+I(XrayRatio*XrayRatio), data = data_rls, nvmax = 13,
                     method = "seqrep")
summary(m.sr)
## Subset selection object
## Call: regsubsets.formula(LengthStay ~ . + I(InfectRisk * InfectRisk) + 
##     I(XrayRatio * XrayRatio), data = data_rls, nvmax = 13, method = "seqrep")
## 14 Variables  (and intercept)
##                            Forced in Forced out
## InfectRisk                     FALSE      FALSE
## Facilities                     FALSE      FALSE
## XrayRatio                      FALSE      FALSE
## Age                            FALSE      FALSE
## CultRatio                      FALSE      FALSE
## NBeds                          FALSE      FALSE
## MedSchool1                     FALSE      FALSE
## Region2                        FALSE      FALSE
## Region3                        FALSE      FALSE
## Region4                        FALSE      FALSE
## DailyCensus                    FALSE      FALSE
## NNurses                        FALSE      FALSE
## I(InfectRisk * InfectRisk)     FALSE      FALSE
## I(XrayRatio * XrayRatio)       FALSE      FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: 'sequential replacement'
##           InfectRisk Facilities XrayRatio Age CultRatio NBeds MedSchool1
## 1  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 2  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 3  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 4  ( 1 )  " "        " "        " "       " " " "       "*"   " "       
## 5  ( 1 )  " "        " "        " "       "*" " "       "*"   " "       
## 6  ( 1 )  " "        " "        " "       "*" " "       "*"   " "       
## 7  ( 1 )  "*"        "*"        "*"       "*" "*"       "*"   "*"       
## 8  ( 1 )  " "        " "        " "       "*" " "       " "   " "       
## 9  ( 1 )  "*"        "*"        "*"       "*" "*"       "*"   "*"       
## 10  ( 1 ) "*"        "*"        "*"       "*" "*"       "*"   "*"       
## 11  ( 1 ) "*"        " "        "*"       "*" " "       "*"   " "       
## 12  ( 1 ) "*"        "*"        "*"       "*" "*"       "*"   "*"       
## 13  ( 1 ) "*"        "*"        "*"       "*" "*"       "*"   "*"       
##           Region2 Region3 Region4 DailyCensus NNurses
## 1  ( 1 )  " "     " "     " "     " "         " "    
## 2  ( 1 )  " "     " "     "*"     " "         " "    
## 3  ( 1 )  " "     " "     "*"     "*"         " "    
## 4  ( 1 )  " "     " "     "*"     "*"         " "    
## 5  ( 1 )  " "     " "     "*"     "*"         " "    
## 6  ( 1 )  " "     " "     "*"     "*"         " "    
## 7  ( 1 )  " "     " "     " "     " "         " "    
## 8  ( 1 )  "*"     "*"     "*"     "*"         "*"    
## 9  ( 1 )  "*"     "*"     " "     " "         " "    
## 10  ( 1 ) "*"     "*"     "*"     " "         " "    
## 11  ( 1 ) "*"     "*"     "*"     "*"         "*"    
## 12  ( 1 ) "*"     "*"     "*"     "*"         "*"    
## 13  ( 1 ) "*"     "*"     "*"     "*"         "*"    
##           I(InfectRisk * InfectRisk) I(XrayRatio * XrayRatio)
## 1  ( 1 )  "*"                        " "                     
## 2  ( 1 )  "*"                        " "                     
## 3  ( 1 )  "*"                        " "                     
## 4  ( 1 )  "*"                        " "                     
## 5  ( 1 )  "*"                        " "                     
## 6  ( 1 )  "*"                        "*"                     
## 7  ( 1 )  " "                        " "                     
## 8  ( 1 )  "*"                        "*"                     
## 9  ( 1 )  " "                        " "                     
## 10  ( 1 ) " "                        " "                     
## 11  ( 1 ) "*"                        "*"                     
## 12  ( 1 ) " "                        " "                     
## 13  ( 1 ) "*"                        " "

Consideramos a inclusão de 1 até o total das 11 variáveis explicativas disponíveis. Em cada um desses casos, o resumo identifica as variáveis incluídas no “melhor” sub-modelo de cada tamanho. No entanto, não temos uma medida da qualidade dos modelos. Precisamos comparar o desempenho dos diferentes modelos obtidos. Para isso, vamos utilizar validação cruzada k-fold (leave-k-out cross-validation):

# Carrega pacote necessário
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
set.seed(111)

train.set <- trainControl(method = "cv", number = 10)
# Calibração do modelo
m.step <- train(LengthStay ~. +I(InfectRisk*InfectRisk)+I(XrayRatio*XrayRatio), data = data_rls, 
                method = "leapBackward", 
                tuneGrid = data.frame(nvmax = 1:13),
                trControl = train.set)
(results.sr <- m.step$results)

Podemos, ainda, selecionar o modelo com melhor ajuste (segundo as métricas consideradas) da seguinte maneira:

m.step$bestTune

Poderíamos, também, escolher o modelo mais reduzido que apresentou melhora comparado aos outros (maior redução no RMSE e MAE e maior aumento no R2). Para visualizar melhor, faremos uma série de plots

# Visualização
par(mfrow=c(3,1))
plot(RMSE ~ nvmax, data = results.sr, type = "l", col = 2, lwd = 2)
plot(Rsquared ~ nvmax, data = results.sr, type = "l", col = 2, lwd = 2)
plot(MAE ~ nvmax, data = results.sr, type = "l", col = 2, lwd = 2)

Nota-se que a adição de 8 variáveis possui o maior impacto na queda de RMSE e MAE e no aumento de Raj^2. As variáveis explicativas são Age, Region, DailyCensus, NNurses, InfectRisk^2 e XrayRatio^2. Vamos considerar como final, o modelo com 8 variáveis explicativas(incluindo as region dummy), pois nota-se que a adição de mais variáveis além de 8 não altera a qualidade do modelo, inclusive a adição de mais de 10 variáveis piora o modelo. Optou-se pelo “melhor” modelo com menor custo. Os coeficientes estimados para o modelo final (8) podem ser recuperados da seguinte maneira:

coef(m.step$finalModel, 8)
##                (Intercept)                        Age 
##               4.742895e+00               6.505347e-02 
##                    Region2                    Region3 
##              -9.176004e-01              -1.199390e+00 
##                    Region4                DailyCensus 
##              -1.931493e+00               9.587186e-03 
##                    NNurses I(InfectRisk * InfectRisk) 
##              -7.005950e-03               5.383361e-02 
##   I(XrayRatio * XrayRatio) 
##               8.965803e-05

PARTE 2

Faremos a all-possible-regressions para identificar um conjunto de variáveis explicativas para o modelo de regressão para a resposta LengthStay. Faremos o modelo completo m.all para aplicar o método AIC.

m.aic <- step(m.all)
## Start:  AIC=57.61
## LengthStay ~ InfectRisk + Facilities + XrayRatio + Age + CultRatio + 
##     NBeds + MedSchool + Region + DailyCensus + NNurses + I(InfectRisk * 
##     InfectRisk) + I(XrayRatio * XrayRatio)
## 
##                              Df Sum of Sq    RSS    AIC
## - MedSchool                   1    0.0975 144.37 55.682
## - CultRatio                   1    0.0994 144.37 55.684
## - InfectRisk                  1    0.3159 144.59 55.853
## - Facilities                  1    0.3825 144.65 55.905
## - XrayRatio                   1    0.6628 144.93 56.124
## - I(XrayRatio * XrayRatio)    1    1.4817 145.75 56.761
## <none>                                    144.27 57.606
## - I(InfectRisk * InfectRisk)  1    2.8734 147.14 57.834
## - NBeds                       1    3.9050 148.18 58.624
## - Age                         1    7.5964 151.87 61.405
## - NNurses                     1    9.3606 153.63 62.710
## - Region                      3   25.6943 169.97 70.127
## - DailyCensus                 1   19.9858 164.26 70.266
## 
## Step:  AIC=55.68
## LengthStay ~ InfectRisk + Facilities + XrayRatio + Age + CultRatio + 
##     NBeds + Region + DailyCensus + NNurses + I(InfectRisk * InfectRisk) + 
##     I(XrayRatio * XrayRatio)
## 
##                              Df Sum of Sq    RSS    AIC
## - CultRatio                   1    0.1505 144.52 53.800
## - Facilities                  1    0.3357 144.70 53.945
## - InfectRisk                  1    0.4154 144.78 54.007
## - XrayRatio                   1    0.6524 145.02 54.192
## - I(XrayRatio * XrayRatio)    1    1.4713 145.84 54.828
## <none>                                    144.37 55.682
## - I(InfectRisk * InfectRisk)  1    3.1538 147.52 56.124
## - NBeds                       1    4.3253 148.69 57.018
## - Age                         1    7.5064 151.87 59.410
## - NNurses                     1    9.2680 153.64 60.713
## - Region                      3   25.7807 170.15 68.249
## - DailyCensus                 1   22.8453 167.21 70.282
## 
## Step:  AIC=53.8
## LengthStay ~ InfectRisk + Facilities + XrayRatio + Age + NBeds + 
##     Region + DailyCensus + NNurses + I(InfectRisk * InfectRisk) + 
##     I(XrayRatio * XrayRatio)
## 
##                              Df Sum of Sq    RSS    AIC
## - Facilities                  1    0.3484 144.87 52.072
## - InfectRisk                  1    0.4791 145.00 52.174
## - XrayRatio                   1    0.5519 145.07 52.231
## - I(XrayRatio * XrayRatio)    1    1.3491 145.87 52.850
## <none>                                    144.52 53.800
## - I(InfectRisk * InfectRisk)  1    3.6952 148.21 54.653
## - NBeds                       1    4.2092 148.73 55.044
## - Age                         1    7.6658 152.18 57.640
## - NNurses                     1    9.1269 153.65 58.720
## - Region                      3   28.8754 173.39 68.384
## - DailyCensus                 1   22.8907 167.41 68.415
## 
## Step:  AIC=52.07
## LengthStay ~ InfectRisk + XrayRatio + Age + NBeds + Region + 
##     DailyCensus + NNurses + I(InfectRisk * InfectRisk) + I(XrayRatio * 
##     XrayRatio)
## 
##                              Df Sum of Sq    RSS    AIC
## - XrayRatio                   1    0.4958 145.36 50.458
## - InfectRisk                  1    0.7351 145.60 50.644
## - I(XrayRatio * XrayRatio)    1    1.2716 146.14 51.060
## <none>                                    144.87 52.072
## - I(InfectRisk * InfectRisk)  1    4.3858 149.25 53.442
## - NBeds                       1    5.2233 150.09 54.075
## - Age                         1    7.5085 152.38 55.782
## - NNurses                     1    9.8972 154.76 57.540
## - Region                      3   28.8744 173.74 66.610
## - DailyCensus                 1   23.8051 168.67 67.264
## 
## Step:  AIC=50.46
## LengthStay ~ InfectRisk + Age + NBeds + Region + DailyCensus + 
##     NNurses + I(InfectRisk * InfectRisk) + I(XrayRatio * XrayRatio)
## 
##                              Df Sum of Sq    RSS    AIC
## - InfectRisk                  1    1.5500 146.91 49.657
## <none>                                    145.36 50.458
## - NBeds                       1    5.1919 150.55 52.424
## - I(InfectRisk * InfectRisk)  1    6.5549 151.92 53.442
## - I(XrayRatio * XrayRatio)    1    7.5029 152.87 54.145
## - Age                         1    7.8623 153.22 54.411
## - NNurses                     1    9.4276 154.79 55.559
## - Region                      3   28.7766 174.14 64.869
## - DailyCensus                 1   23.6348 169.00 65.482
## 
## Step:  AIC=49.66
## LengthStay ~ Age + NBeds + Region + DailyCensus + NNurses + I(InfectRisk * 
##     InfectRisk) + I(XrayRatio * XrayRatio)
## 
##                              Df Sum of Sq    RSS    AIC
## <none>                                    146.91 49.657
## - NBeds                       1    4.7261 151.64 51.235
## - I(XrayRatio * XrayRatio)    1    6.8032 153.72 52.772
## - Age                         1    9.4577 156.37 54.707
## - NNurses                     1   11.1109 158.02 55.895
## - DailyCensus                 1   22.8907 169.80 64.019
## - Region                      3   30.2693 177.18 64.826
## - I(InfectRisk * InfectRisk)  1   29.3250 176.24 68.222

O teste indica que o melhor modelo deve ser feito usando as seguintes variáveis: NBeds, I(XrayRatio * XrayRatio), Age, NNuurses, DailyCensus, Region e I(InfectRisk * InfectRisk). Seguiremos usando o critério para decisão. Poderíamos realizar o mesmo teste da seguinte maneira:

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
m.aic2 <- stepAIC(m.all, direction = "both", trace = FALSE)
summary(m.aic2)
## 
## Call:
## lm(formula = LengthStay ~ Age + NBeds + Region + DailyCensus + 
##     NNurses + I(InfectRisk * InfectRisk) + I(XrayRatio * XrayRatio), 
##     data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3187 -0.7407 -0.0855  0.5983  5.9795 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4.626e+00  1.485e+00   3.115 0.002383 ** 
## Age                         6.720e-02  2.610e-02   2.575 0.011443 *  
## NBeds                      -6.022e-03  3.308e-03  -1.820 0.071620 .  
## Region2                    -7.798e-01  3.297e-01  -2.365 0.019885 *  
## Region3                    -1.091e+00  3.261e-01  -3.346 0.001147 ** 
## Region4                    -1.770e+00  4.029e-01  -4.393 2.72e-05 ***
## DailyCensus                 1.611e-02  4.021e-03   4.006 0.000117 ***
## NNurses                    -5.836e-03  2.091e-03  -2.791 0.006262 ** 
## I(InfectRisk * InfectRisk)  5.249e-02  1.158e-02   4.534 1.57e-05 ***
## I(XrayRatio * XrayRatio)    8.971e-05  4.107e-05   2.184 0.031234 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.194 on 103 degrees of freedom
## Multiple R-squared:  0.641,  Adjusted R-squared:  0.6096 
## F-statistic: 20.43 on 9 and 103 DF,  p-value: < 2.2e-16

Note que o resultado apresentou as mesmas 7 variáveis(incluindo as dummy region).

Cp de Mallows e R2aj

library(leaps)

data_rls_maisq <- cbind(data_rls, IR2 = InfectRisk*InfectRisk, XR2 = XrayRatio*XrayRatio)
data_rls_maisq <- data_rls_maisq[-c(47,112), ]
attach(data_rls_maisq)
## The following objects are masked from scenic:
## 
##     Age, CultRatio, DailyCensus, Facilities, InfectRisk, LengthStay,
##     MedSchool, NBeds, NNurses, Region, XrayRatio
m.crit <- regsubsets(LengthStay ~., data = na.omit(data_rls_maisq), nvmax = 13,
                     method = "exhaustive")
(results.crit <- summary(m.crit))
## Subset selection object
## Call: regsubsets.formula(LengthStay ~ ., data = na.omit(data_rls_maisq), 
##     nvmax = 13, method = "exhaustive")
## 14 Variables  (and intercept)
##             Forced in Forced out
## InfectRisk      FALSE      FALSE
## Facilities      FALSE      FALSE
## XrayRatio       FALSE      FALSE
## Age             FALSE      FALSE
## CultRatio       FALSE      FALSE
## NBeds           FALSE      FALSE
## MedSchool1      FALSE      FALSE
## Region2         FALSE      FALSE
## Region3         FALSE      FALSE
## Region4         FALSE      FALSE
## DailyCensus     FALSE      FALSE
## NNurses         FALSE      FALSE
## IR2             FALSE      FALSE
## XR2             FALSE      FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: exhaustive
##           InfectRisk Facilities XrayRatio Age CultRatio NBeds MedSchool1
## 1  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 2  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 3  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 4  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 5  ( 1 )  " "        " "        " "       " " " "       " "   " "       
## 6  ( 1 )  " "        " "        " "       "*" " "       " "   " "       
## 7  ( 1 )  " "        " "        " "       "*" " "       " "   " "       
## 8  ( 1 )  " "        " "        " "       "*" " "       " "   " "       
## 9  ( 1 )  " "        " "        " "       "*" " "       " "   "*"       
## 10  ( 1 ) " "        " "        "*"       "*" " "       " "   "*"       
## 11  ( 1 ) " "        "*"        "*"       "*" " "       " "   "*"       
## 12  ( 1 ) " "        "*"        "*"       "*" " "       "*"   "*"       
## 13  ( 1 ) " "        "*"        "*"       "*" "*"       "*"   "*"       
##           Region2 Region3 Region4 DailyCensus NNurses IR2 XR2
## 1  ( 1 )  " "     " "     " "     " "         " "     "*" " "
## 2  ( 1 )  " "     " "     "*"     " "         " "     "*" " "
## 3  ( 1 )  " "     " "     "*"     "*"         " "     "*" " "
## 4  ( 1 )  " "     " "     "*"     "*"         " "     "*" "*"
## 5  ( 1 )  "*"     "*"     "*"     "*"         " "     "*" " "
## 6  ( 1 )  "*"     "*"     "*"     "*"         " "     "*" " "
## 7  ( 1 )  "*"     "*"     "*"     "*"         " "     "*" "*"
## 8  ( 1 )  "*"     "*"     "*"     "*"         "*"     "*" "*"
## 9  ( 1 )  "*"     "*"     "*"     "*"         "*"     "*" "*"
## 10  ( 1 ) "*"     "*"     "*"     "*"         "*"     "*" "*"
## 11  ( 1 ) "*"     "*"     "*"     "*"         "*"     "*" "*"
## 12  ( 1 ) "*"     "*"     "*"     "*"         "*"     "*" "*"
## 13  ( 1 ) "*"     "*"     "*"     "*"         "*"     "*" "*"
n.var <- 1:13 # numero de var. explicativas
plot(n.var + 1, results.crit$cp, xlab = "no. parâmetros", ylab = "Cp de Mallows")
abline(0,1, col = "red")

Queremos modelos que estejam abaixo da linha \(C_p=p+1\), em que p é o no. de variáveis explicativas no modelo. Neste caso, o melhor modelo é o que inclui 9 variáveis explicativas, contando as Regions dummy. As variáveis são Age, NBeds, Region, DailyCensus, NNurses, IR2 e XR2.

Vejamos os resultados utilizando como critério de avaliação o coeficiente de determinação ajustado:

n.var <- 1:13 # numero de var. explicativas
plot(n.var + 1, results.crit$adjr2, 
     xlab = "no. parâmetros", 
     ylab = "R2 ajustado", type="b")

Percebe-se que nesse caso, o melhor modelo também inclui 9 variáveis(incluindo as Regions dummy). Nota-se que consistentemente as variáveis Region, DailyCensus, NNurses, IR2(a qual pelo principio da hierarquia requer IR no modelo) foram selecionadas por todos os métodos. Além disso, tem-se que Age, XR2 figuram nos outros 3 métodos sugeridos. Portanto, a fim de analisar a interação dessas variáveis (Age e XR2), será construído um modelo sugerido pelo Stepwise que inclui as variáveis que aparecem em todos métodos além de Age e XR2.

Partiremos agora para a construção e análise do modelo com essas variáveis selecionadas para a resposta LengthStay:

Vamos construir um modelo de 2a. ordem (resposta quadrática) para a variável LengthStay como função de variáveis explicativas, mas com atenção especial para IR e XR que tem termos quadráticos.

Vamos verificar a correlação entre a variável explicativa IR e IR2 e XR e XR2. Caso a correlação seja elevada, é desejável centralizar a variável explicativa quantitativa, fazendo InfectRisk.c = InfectRisk - mean(InfectRisk) e XrayRatio.c = XrayRatio - mean(XrayRatio).

cor(InfectRisk, IR2)
## [1] 0.9745198
cor(XrayRatio, XR2)
## [1] 0.9880858

A correlação foi elevada, portanto será feita a centralização em torno da média.

InfectRisk.c <- InfectRisk - mean(InfectRisk)
XrayRatio.c <- XrayRatio - mean(XrayRatio)
data_maisq_cent <- data.frame(LengthStay = LengthStay, InfectRisk.c = InfectRisk.c, Facilities = Facilities, XrayRatio.c = XrayRatio.c, Age = Age, CultRatio =  CultRatio, NBeds =  NBeds, MedSchool = MedSchool, Region = Region, DailyCensus = DailyCensus, NNurses = NNurses, IR2.c = InfectRisk.c*InfectRisk.c, XR2.c = XrayRatio.c*XrayRatio.c )

attach(data_maisq_cent)
## The following objects are masked _by_ .GlobalEnv:
## 
##     InfectRisk.c, XrayRatio.c
## The following objects are masked from data_rls_maisq:
## 
##     Age, CultRatio, DailyCensus, Facilities, LengthStay, MedSchool,
##     NBeds, NNurses, Region
## The following objects are masked from scenic:
## 
##     Age, CultRatio, DailyCensus, Facilities, LengthStay, MedSchool,
##     NBeds, NNurses, Region
modelo_rlm <- lm(LengthStay ~ Age + DailyCensus + NNurses + IR2.c + XR2.c + InfectRisk.c + XrayRatio.c, data  = data_maisq_cent)

summary(modelo_rlm)
## 
## Call:
## lm(formula = LengthStay ~ Age + DailyCensus + NNurses + IR2.c + 
##     XR2.c + InfectRisk.c + XrayRatio.c, data = data_maisq_cent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.24279 -0.84708  0.04188  0.52402  2.88169 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.3811807  1.3158616   4.849 4.40e-06 ***
## Age           0.0438328  0.0245323   1.787  0.07692 .  
## DailyCensus   0.0058456  0.0019908   2.936  0.00410 ** 
## NNurses      -0.0030747  0.0021255  -1.447  0.15106    
## IR2.c         0.0884198  0.0451966   1.956  0.05313 .  
## XR2.c         0.0001694  0.0002384   0.711  0.47885    
## InfectRisk.c  0.4023512  0.0974693   4.128 7.45e-05 ***
## XrayRatio.c   0.0179353  0.0062391   2.875  0.00491 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.115 on 103 degrees of freedom
## Multiple R-squared:  0.4644, Adjusted R-squared:  0.428 
## F-statistic: 12.76 on 7 and 103 DF,  p-value: 1.051e-11

A abordagem utilizada na construção do modelo consistirá na criação de um modelo inicial incluindo apenas as variáveis quantitativas, já considerando a centralização das variáveis InfectRisk e XrayRatio. Em seguida, após a análise deste modelo e remoção de eventuais variáveis, serão incluídos e analisados os termos de interação entre as variáveis quantitativas e a variável categórica Region.

plot(modelo_rlm)

Devido à baixa significância de XR2 e de NNurses no modelo feito, e para simplificação, XR2 e NNurses será removido

modelo_rlm2 <- lm(LengthStay ~ Age + DailyCensus + IR2.c + InfectRisk.c + XrayRatio.c, data  = data_maisq_cent)
summary(modelo_rlm2)
## 
## Call:
## lm(formula = LengthStay ~ Age + DailyCensus + IR2.c + InfectRisk.c + 
##     XrayRatio.c, data = data_maisq_cent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.22656 -0.83740  0.04597  0.52906  2.93223 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.2779409  1.3142453   4.777 5.80e-06 ***
## Age          0.0452076  0.0244766   1.847 0.067567 .  
## DailyCensus  0.0033236  0.0008157   4.074 8.98e-05 ***
## IR2.c        0.1075181  0.0422021   2.548 0.012292 *  
## InfectRisk.c 0.3826969  0.0967448   3.956 0.000139 ***
## XrayRatio.c  0.0180090  0.0062382   2.887 0.004724 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.117 on 105 degrees of freedom
## Multiple R-squared:  0.4524, Adjusted R-squared:  0.4263 
## F-statistic: 17.35 on 5 and 105 DF,  p-value: 1.735e-12
plot(modelo_rlm2)

Em seguida analisaremos se há indícios de multicolinearidade entre as variáveis explicativas.

data_multicol <- data.frame(LengthStay = LengthStay, Age = Age, DailyCensus = DailyCensus, IR2.c = IR2.c, InfectRisk.c = InfectRisk.c, XrayRatio.c = XrayRatio.c)

attach(data_multicol)
## The following objects are masked _by_ .GlobalEnv:
## 
##     InfectRisk.c, XrayRatio.c
## The following objects are masked from data_maisq_cent:
## 
##     Age, DailyCensus, InfectRisk.c, IR2.c, LengthStay, XrayRatio.c
## The following objects are masked from data_rls_maisq:
## 
##     Age, DailyCensus, LengthStay
## The following objects are masked from scenic:
## 
##     Age, DailyCensus, LengthStay
modelo_rlm_multicol <- lm(LengthStay ~ Age + DailyCensus + IR2.c + InfectRisk.c + XrayRatio.c, data  = data_multicol)
summary(modelo_rlm_multicol)
## 
## Call:
## lm(formula = LengthStay ~ Age + DailyCensus + IR2.c + InfectRisk.c + 
##     XrayRatio.c, data = data_multicol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.22656 -0.83740  0.04597  0.52906  2.93223 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.2779409  1.3142453   4.777 5.80e-06 ***
## Age          0.0452076  0.0244766   1.847 0.067567 .  
## DailyCensus  0.0033236  0.0008157   4.074 8.98e-05 ***
## IR2.c        0.1075181  0.0422021   2.548 0.012292 *  
## InfectRisk.c 0.3826969  0.0967448   3.956 0.000139 ***
## XrayRatio.c  0.0180090  0.0062382   2.887 0.004724 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.117 on 105 degrees of freedom
## Multiple R-squared:  0.4524, Adjusted R-squared:  0.4263 
## F-statistic: 17.35 on 5 and 105 DF,  p-value: 1.735e-12
plot(modelo_rlm_multicol)

cor(data_multicol)
##              LengthStay         Age DailyCensus       IR2.c InfectRisk.c
## LengthStay    1.0000000  0.12220505  0.41582956  0.12007536   0.54946732
## Age           0.1222051  1.00000000 -0.09364990  0.19343483  -0.02842443
## DailyCensus   0.4158296 -0.09364990  1.00000000 -0.16256630   0.36526980
## IR2.c         0.1200754  0.19343483 -0.16256630  1.00000000  -0.06803173
## InfectRisk.c  0.5494673 -0.02842443  0.36526980 -0.06803173   1.00000000
## XrayRatio.c   0.3763310 -0.04563341  0.03935332 -0.08637082   0.43755783
##              XrayRatio.c
## LengthStay    0.37633102
## Age          -0.04563341
## DailyCensus   0.03935332
## IR2.c        -0.08637082
## InfectRisk.c  0.43755783
## XrayRatio.c   1.00000000
corrplot(cor(data_multicol), method = "circle", order = "hclust", sig.level = 0.01, insig = "blank")

vif(modelo_rlm_multicol)
##          Age  DailyCensus        IR2.c InfectRisk.c  XrayRatio.c 
##     1.044571     1.214626     1.068809     1.457343     1.274040

Uma observação inicial dos gráficos de dispersão obtidos sugere que não há aparente relação linear entre as variáveis explicativas. Em seguida, a matriz de correlação corrobora essa hipótese, uma vez que não apresentou nenhum valor próximo de 1 (o valor que mais se aproxima de 1 foi 0.43755783 para InfectRisk.c versus XrayRatio.c). Em seguida, fez-se o VIF para o modelo criado. Para esse teste, temos que valores acima de 5 indicam sérios efeitos de multicolinearidade, sendo que o coeficiente de regressão para a variável não está sendo estimado de maneira apropriada. Para o modelo obtido foram encontrados valores próximos de 1, indicando que não há aumento relevante da variância dos coeficientes da regressão devido à presença de multicolinearidade. Dessa maneira, os três resultados obtidos indicam que não há efeitos sérios de multicolinearidade presentes no modelo construido, de modo que não serão considerados os termos de interação entre as variáveis quantitativas.

Será adicionada a variável categórica Region ao modelo.

modelo.quali <- lm(LengthStay ~ Age + Region + DailyCensus + IR2.c + InfectRisk.c + XrayRatio.c)
summary(modelo.quali)
## 
## Call:
## lm(formula = LengthStay ~ Age + Region + DailyCensus + IR2.c + 
##     InfectRisk.c + XrayRatio.c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.13526 -0.68976 -0.02954  0.51331  3.01634 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.2150526  1.2118351   5.954 3.74e-08 ***
## Age           0.0441287  0.0220349   2.003 0.047866 *  
## Region2      -0.5671194  0.2686772  -2.111 0.037236 *  
## Region3      -0.8520622  0.2675502  -3.185 0.001922 ** 
## Region4      -1.8374216  0.3271966  -5.616 1.70e-07 ***
## DailyCensus   0.0027013  0.0007445   3.628 0.000448 ***
## IR2.c         0.0777306  0.0382318   2.033 0.044636 *  
## InfectRisk.c  0.4034252  0.0873062   4.621 1.12e-05 ***
## XrayRatio.c   0.0103694  0.0056919   1.822 0.071418 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9865 on 102 degrees of freedom
## Multiple R-squared:  0.5851, Adjusted R-squared:  0.5525 
## F-statistic: 17.98 on 8 and 102 DF,  p-value: < 2.2e-16
plot(modelo.quali)

Pelo p-valor ter diminuído e o valor de R² ajustado ter aumentado em relação ao modelo sem a variável categórica Region, então o acréscimo de tal variável é relevante ao nosso modelo.

Agora, analisaremos as interação entre os termos da variável categórica e as variáveis quantitativas.

modelo_rlm_cat <- lm(LengthStay ~ Age + Region + DailyCensus + IR2.c + InfectRisk.c + XrayRatio.c + Region*Age + Region*DailyCensus + Region*IR2.c + (Region*InfectRisk.c) + (Region*XrayRatio.c), data  = data_maisq_cent)
summary(modelo_rlm_cat)
## 
## Call:
## lm(formula = LengthStay ~ Age + Region + DailyCensus + IR2.c + 
##     InfectRisk.c + XrayRatio.c + Region * Age + Region * DailyCensus + 
##     Region * IR2.c + (Region * InfectRisk.c) + (Region * XrayRatio.c), 
##     data = data_maisq_cent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.36469 -0.56340 -0.00599  0.53684  3.05013 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           3.3533086  3.0007688   1.117  0.26686   
## Age                   0.1175777  0.0566188   2.077  0.04078 * 
## Region2               4.5712149  3.7369374   1.223  0.22454   
## Region3               2.5287898  3.6849185   0.686  0.49438   
## Region4               4.2441957  4.1383413   1.026  0.30793   
## DailyCensus           0.0027308  0.0018616   1.467  0.14601   
## IR2.c                -0.0582324  0.1106111  -0.526  0.59991   
## InfectRisk.c          0.6478860  0.2322439   2.790  0.00648 **
## XrayRatio.c           0.0124189  0.0122559   1.013  0.31373   
## Age:Region2          -0.1035561  0.0701189  -1.477  0.14332   
## Age:Region3          -0.0496499  0.0693362  -0.716  0.47586   
## Age:Region4          -0.1260798  0.0786987  -1.602  0.11277   
## Region2:DailyCensus   0.0008952  0.0022425   0.399  0.69071   
## Region3:DailyCensus  -0.0031308  0.0023292  -1.344  0.18240   
## Region4:DailyCensus   0.0024598  0.0028433   0.865  0.38934   
## Region2:IR2.c         0.1946734  0.1282904   1.517  0.13278   
## Region3:IR2.c         0.0493336  0.1299974   0.379  0.70524   
## Region4:IR2.c         0.5191157  0.3596158   1.444  0.15246   
## Region2:InfectRisk.c -0.3067356  0.2781482  -1.103  0.27316   
## Region3:InfectRisk.c -0.0971259  0.2777348  -0.350  0.72740   
## Region4:InfectRisk.c -0.5165104  0.4002669  -1.290  0.20033   
## Region2:XrayRatio.c  -0.0059920  0.0162800  -0.368  0.71372   
## Region3:XrayRatio.c  -0.0092088  0.0184164  -0.500  0.61831   
## Region4:XrayRatio.c  -0.0049160  0.0193012  -0.255  0.79955   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9746 on 87 degrees of freedom
## Multiple R-squared:  0.6545, Adjusted R-squared:  0.5632 
## F-statistic: 7.167 on 23 and 87 DF,  p-value: 5.372e-12
anova(modelo_rlm_cat)

Olhando a coluna do valor F, podemos perceber que as interações de Region com as variáveis quantitativas possuem um valor F pequeno, indicando que podemos descartá-las. Isso também pode ser verificado pela comparação entre os dois modelos a seguir.

anova(modelo.quali, modelo_rlm_cat, test="F")
modelo.quali
## 
## Call:
## lm(formula = LengthStay ~ Age + Region + DailyCensus + IR2.c + 
##     InfectRisk.c + XrayRatio.c)
## 
## Coefficients:
##  (Intercept)           Age       Region2       Region3       Region4  
##     7.215053      0.044129     -0.567119     -0.852062     -1.837422  
##  DailyCensus         IR2.c  InfectRisk.c   XrayRatio.c  
##     0.002701      0.077731      0.403425      0.010369

Analisando os resíduos do modelo, obtemos:

summary(modelo.quali)
## 
## Call:
## lm(formula = LengthStay ~ Age + Region + DailyCensus + IR2.c + 
##     InfectRisk.c + XrayRatio.c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.13526 -0.68976 -0.02954  0.51331  3.01634 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.2150526  1.2118351   5.954 3.74e-08 ***
## Age           0.0441287  0.0220349   2.003 0.047866 *  
## Region2      -0.5671194  0.2686772  -2.111 0.037236 *  
## Region3      -0.8520622  0.2675502  -3.185 0.001922 ** 
## Region4      -1.8374216  0.3271966  -5.616 1.70e-07 ***
## DailyCensus   0.0027013  0.0007445   3.628 0.000448 ***
## IR2.c         0.0777306  0.0382318   2.033 0.044636 *  
## InfectRisk.c  0.4034252  0.0873062   4.621 1.12e-05 ***
## XrayRatio.c   0.0103694  0.0056919   1.822 0.071418 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9865 on 102 degrees of freedom
## Multiple R-squared:  0.5851, Adjusted R-squared:  0.5525 
## F-statistic: 17.98 on 8 and 102 DF,  p-value: < 2.2e-16

Aqui podemos verificar que o p-valor do modelo é baixo e o R² ajustado, como vimos no passo a passo, apenas aumentou com a adição das variáveis explciativas.

plot(modelo.quali)

O gráfico de resíduos possui uma leve inclinação em relação à horizontal. Assim, realizaremos o teste de Durbin-Watson.

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
dwtest(modelo.quali, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  modelo.quali
## DW = 1.9339, p-value = 0.7152
## alternative hypothesis: true autocorrelation is not 0

Portanto, devido aos valores-p encontrados, aceita-se a hipótese nula, ou seja, sugere-se que os erros não são correlacionados.

E realizaremos também o teste de Shapiro-Wilk.

resid_rls_pad_ir <- rstandard(modelo.quali)
shapiro.test(resid_rls_pad_ir)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_rls_pad_ir
## W = 0.98253, p-value = 0.1556

O valor-p do teste nos leva a aceitar a hipótese nula e concluir que os erros são normalmente distribuídos. Tais resultados indicam que não é necessário realizar transformações e nem realizar outros tipos de regressão especial. A leve inclinação no gráfico de resíduos provavelmente está relacionada ao fato de que os termos de interação entre variáveis foram considerados desprezíveis, contudo ainda existe uma pequena influência que pode ser responsável.

Conclusão

A base de dados foi inicialmente submetida a uma análise exploratória preliminar e em seguida aplicados os métodos de seleção de variáveis para escolha do modelo, sendo escolhido o modelo proposto pelo método Stepwise Regression, pois apresentou resultados coerentes com os demais métodos. Vale ressaltar que foi incluída um termo de segunda ordem nessa etapa do processo, consistindo de InfectRisk².

Em seguida, foi construído um modelo com as variáveis identificadas, partindo-se de um modelo inicial com apenas variáveis quantitativas, no qual foi discutido e descartado a possível existência de multicolinearidade entre elas. Após isso, adicionou-se a variável categórica e analisou-se a sua influência no modelo, bem como termos de interação entre ela e as quantitativas. Após descartado os coeficientes não significativos, obtemos um modelo final para explicar a variável LenghtStay que consiste das variáveis explicativas Age, Region, DailyCensus, IR2.c, InfectRisk.c, XrayRatio.c, dado por:

\(E[LenghtStay]=β_0+β_1*Age+β_2*Region+β_3*DailyCensus+β_4*IR2.c+β_5*InfectRisk.c+β_6*XrayRatio.c\)

Dado que LenghtStay é o período de internação médio de todos os pacientes no hospital (em dias), podemos concluir que tal variável pode ser explicado por Age (que é a idade média dos pacientes), Region (região geográfica), DailyCensus (número médio de pacientes no hospital por dia, durante o período avaliado), XrayRatio (razão do número de raios-X realizados pelo número de pacientes sem sintomas de pneumonia), InfectRisk (risco de infecção, calculado como a probabilidade média estimada de contrair infecção no hospital e InfectRisk². Note que esse modelo é linear para todas as variáveis exceto InfectRisk. Algumas considerações importantes podem ser feitos com isso: o modelo obtido possui grau polinomial de baixa ordem, de modo que não há ruídos aletórios sendo modelados (fenômeno de Rudge). Além disso é importante ressaltar a característica da interação, uma vez que o efeito da variável InfectRisk na resposta depende dessa variável.

\(\partial{E[Y]}/\partial{(InfectRisk.c)} = β_5 + 2β_4*InfectRisk.c\)

Ademais, para o modelo polinomial ser utilizado para inferências ou testes de hipóteses é muito importante que seja feito dentro do escopo, de modo que esse modelo não permita extrapolações. Com isso, nota-se que essa variável é importante na previsão do período médio de internação de um paciente LengthStay. Assim, conclui-se que programas de vigilância que conseguem controlar o risco de infecção impactam o período médio de internação.

Os coeficientes R² e o R² ajustado para o modelo não são unitários, de modo que essas variáveis não explicam com totalidade a variável de resposta LenghtStay. Contudo, durante a seleçao de variáveis percebe-se que a inclusão de outras variáveis não apresenta melhoras significativas nesse aspecto. Além disso, o modelo obtido não apresentou indícios de problemas como multicolinearidade, dependências ou heterocedasticidade. DEssa forma, é satisfatório a utilização desse modelo para explicar variações da variável LenghtStay. Por exemplo, podemos utilizar o modelo para verificar a inferência sobre a resposta esperada, desejando estimar o LenghtStay para Age = 50.9, Region = 2, DailyCensus = 147, InfectRisk.c = 0.77837838, InfectRisk² = 6.058729e-01 e XrayRatio = 15.7522523, como mostrado abaixo. Note que o ponto utilizado está dentro do escopo do modelo, coerente com a teoria.

XH <- data.frame(Age = 50.9, Region = '2', DailyCensus = 147, InfectRisk.c = 0.77837838, IR2.c = 6.058729e-01, XrayRatio.c = 15.7522523)
predict.lm(modelo.quali, newdata = XH, interval = "prediction", level = 0.95)
##        fit      lwr      upr
## 1 9.815633 7.817097 11.81417

Concluímos com 95% de confiança que o período de internação médio para Age = 50.9, Region = 2, DailyCensus = 147, InfectRisk.c = 0.77837838, InfectRisk² = 6.058729e-01 e XrayRatio = 15.7522523 deve estar entre 7.817 e 11.814 dias. Dependendo do nível da precisão requerida, o intervalo pode ser adequado para a estimativa. Caso seja necessário um intervalo menor o modelo precisaria de ajustes.